#include "patch.h"
#include <cmath>

inline float w0(float x) { return (1.0f / 6.0f) * (-x*x*x + 3.0f * x*x - 3.0f * x + 1.0f); }
inline float w1(float x) { return (1.0f / 6.0f) * (3.0f * x*x*x - 6.0f * x*x + 4.0f); }
inline float w2(float x) { return (1.0f / 6.0f) * (-3.0f * x*x*x + 3.0f * x*x + 3.0f * x + 1.0f); }
inline float w3(float x) { return (1.0f / 6.0f) * (x*x*x); }

template <class T>
T bicubic(float u, float v, T data[4][4]) {
    float uw[] = { w0(u), w1(u), w2(u), w3(u) };
    float vw[] = { w0(v), w1(v), w2(v), w3(v) };
    T r;

    for (int i = 0; i < 4; ++i) {
        T col;
        for (int j = 0; j < 4; ++j) {
            col += uw[j] * data[i][j];
        }
        r += vw[i] * col;
    }
    return r;
};

void bicubic_patch::split(std::stack<primitive_ptr>& primitives) {
    vertex l[5][5];             //5x5 grid for subdivision step

    //calculate the 9 facepoints by averaging neighboring vertices
    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            l[i*2][j*2] = vertices[i][j] + vertices[i][j+1] + vertices[i+1][j] + vertices[i+1][j+1];
            l[i*2][j*2] /= 4.0f;
        }
    }

    //sum up the edge endpoints around the new edge midpoints
    //average later after adding the sum to the original points
    for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < 3; ++j) { 
            l[i*2+1][j*2] = vertices[i+1][j+0] + vertices[i+1][j+1]; 
            l[j*2][i*2+1] = vertices[j+0][i+1] + vertices[j+1][i+1];
        }
    }

    //calculate the 4 original points
    for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < 2; ++j) {
            //add the original point with it's weight 
            l[i*2+1][j*2+1] = 4.0f * vertices[i+1][j+1];
            vertex m;

            //sum up the 3x3 grid around original point, edge points already had right weight
            for (int ii = 0; ii < 3; ++ii) {
                for (int jj = 0; jj < 3; ++jj) {
                    m += l[ii+i*2][jj+j*2];
                }
            }
            l[i*2+1][j*2+1] = m / 16.0f; //average
        }
    }

    //now calculate edge points by adding neighboring facepoints and averaging
    for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < 3; ++j) {
            l[i*2+1][j*2] += l[i*2+0*2][j*2] + l[i*2+1*2][j*2];
            l[j*2][i*2+1] += l[j*2][i*2+0*2] + l[j*2][i*2+1*2];
            l[i*2+1][j*2] /= 4.0f;
            l[j*2][i*2+1] /= 4.0f;
        }
    }
    
    bicubic_patch* out_patches[4];

    for (int i = 0; i < 4; ++i) {
        out_patches[i] = new bicubic_patch();
    }

    //copy the vertices in output patches
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            out_patches[0]->vertices[i][j] = l[i][j];
            out_patches[1]->vertices[i][j] = l[i+1][j];
            out_patches[2]->vertices[i][j] = l[i][j+1];
            out_patches[3]->vertices[i][j] = l[i+1][j+1];
        }
    }

    for (int i = 0; i < 4; ++i) {
        primitives.push(out_patches[i]);
    }
};


grid<vertex>* bicubic_patch::dice(int uratio, int vratio) {
    grid<vertex>* out = new grid<vertex>(uratio,vratio);

    for (int i = 0; i < uratio; ++i) {
        for (int j = 0; j < vratio; ++j) {
            float u = float(i) / float(uratio - 1);
            float v = float(j) / float(vratio - 1);
            (*out)(i, j) = bicubic(u,v,vertices);
        }
    }
    return out;
}

float triangle_area(const vec2f a, const vec2f b) {
    return 0.5f * fabs(a[0]*b[1] - b[0]*a[1]);
}

inline float fabs(float x) { return (x < 0) ? -x : x; }

bool bicubic_patch::diceable(const matrix4f& view, float shadingrate, int gridsize, int& u, int& v) {
    vec4f tr[4] = { 
        vertices[1][1], vertices[1][2], 
        vertices[2][1], vertices[2][2]
    };
    for (int i = 0; i < 4; ++i) tr[i] = (view * tr[i]).hdivide();
    
    vec2f ab = vec2f(tr[1]) - vec2f(tr[0]);
    vec2f ac = vec2f(tr[2]) - vec2f(tr[0]);
    vec2f db = vec2f(tr[1]) - vec2f(tr[3]);
    vec2f dc = vec2f(tr[2]) - vec2f(tr[3]);

    float el[4] = {
        fabs(ab.norm()), fabs(ac.norm()),
        fabs(db.norm()), fabs(dc.norm())
    };

    for (int i = 0; i < 4; ++i) {
        if (gridsize < (el[i] * 2.0f)) return false;
    }

    float area = triangle_area(ab, ac) + triangle_area(dc, db);

    if (gridsize < area) return false;

    float uu = (el[0] + el[3]) / 2.0f;
    float vv = (el[1] + el[2]) / 2.0f;

    if (gridsize < (uu * vv)) return false;

    u = floor(uu / sqrt(shadingrate));
    v = floor(vv / sqrt(shadingrate));
    
    //grid boundary extrapolation requires center vertices
    if (u < 3) u = 3;
    if (v < 3) v = 3;

    return true;
};
